Impurity-assisted electric control of spin-valley qubits in monolayer MoS₂

The electron spin in confined semiconductor quantum dots [1, 2] (QDs) represents an ideal qubit system for encoding information at the quantum level, and a promising building block for quantum information processing. All-electric manipulation of electron spins in QDs is enabled by spin–orbit interaction in two-dimensional electron systems as well as in nanowires [3–8]. The electronic valley degree of freedom [9], relevant, e.g. in silicon [10, 11], carbon nanostructures, and transition metal dichalcogenides (TMDCs), has also been proposed for quantum information processing purposes. This degree of freedom could be utilized as a qubit on its own [12–14], or together with the electron spin, forming a combined spin-valley qubit [14–18].

The rise of two-dimensional (2D) materials and van der Waals heterostructures [19, 20] promoted 2D semiconducting TMDCs as alternative platforms for electronics, spintronics [21, 22], and valleytronics [23] applications, opening up new opportunities in nanoelectronics and optoelectronics with two-dimensional crystals [24]. An appealing feature of TMDCs is the strong spin–orbit interaction, which is characteristic of both the valence and conduction bands [25, 26], and arises due to the presence of the heavy transition-metal atoms of the material. In particular, the broken inversion symmetry of monolayer (ML) TMDCs gives rise to a strong spin-valley locking, whereby the Bloch states close to the valleys have an out-of-plane spin polarization [27]. Furthermore, the possibility of electrostatically defining QDs in TMDCs such as MoS₂, WS₂ and WSe₂ [28–33] offers new opportunities for spin-based quantum information processing [34–39]. On the one hand, the nuclear-spin-free environment [37, 38], achievable via isotopic purification, is expected to prolong the qubit lifetime compared to III–V materials such as GaAs—in fact, this strategy has already proven to boost decoherence times in diamond and silicon [40]. On the other hand, the spin–orbit interaction present in TMDCs offers the possibility of efficient spin control via electric fields, which has potential advantages over magnetic control.

In this work, we consider a single electron confined in a QD, which is electrostatically defined in the conduction band of a ML TMDC (see figure 1). The two lowest-energy states in this setup form a spin-valley qubit: the strong spin–orbit interaction and the broken inversion symmetry lock the spin to the valley degree of freedom, and thereby the qubit basis states are characterized by opposite spin in opposite valleys, $\renewcommand{\ket}[1]{|#1\rangle} \ket{K \uparrow}$ and $\renewcommand{\ket}[1]{|#1\rangle} \ket{K' \downarrow}$, in analogy with the spin-valley (or 'Kramers') qubit in carbon nanotubes [18]. An interesting feature of the spin-valley qubit is that it is difficult to induce a coupling between the basis states, since that requires a simultaneous flip of the spin and the valley. On the one hand, this is an advantage, since it results in a suppressed qubit relaxation, which might imply a prolonged qubit lifetime. On the other hand, this makes it difficult to control the qubit with resonant excitation.

Zoom In Zoom Out Reset image size

Here, we show that a short-range impurity (e.g. vacancy, substitutional atom, adatom) in a ML-TMDC QD (see figure 1) can couple the basis states of the spin-valley qubit, and thereby allow for resonant qubit control via an ac electric field, in the spirit of electrically driven spin resonance. The two main target quantities we calculate are (1) the in-plane g-factor g_xx, which is made finite by the presence of a short-range impurity, and (2) the Rabi frequency $\Omega_{\rm R}$ characterizing the speed of the electrically induced dynamics of the spin-valley qubit. We calculate how these two quantities depend on the system parameters, in particular their dependence $g_{xx}(x_0, y_0)$ and $\Omega_{\rm R}(x_0, y_0)$ on the relative position of the impurity and the QD centre. Based on recent results [36, 41, 42] on symmetry-forbidden intervalley scattering in ML-TMDCs, we reveal that the spatial patterns $g_{xx}(x_0, y_0)$ and $\Omega_\text{R}(x_0,y_0)$ depend drastically on the type of the impurity (figure 3): qualitatively different spatial patterns are obtained for a defect located at a transition-metal site (Mo-type impurity, e.g. a Mo vacancy in MoS₂), and for one located at a chalcogen site (S-type impurity, e.g. an S vacancy in MoS₂). We discuss the role of symmetry-forbidden intervalley scattering for adatoms as well (see table 1). We also highlight the fact that the impurity-induced coupling between the spin-valley qubit states increases as the spin–orbit splitting decreases. Therefore, we expect that the associated effects are most easily observable in the ML-TMDC with the smallest conduction-band spin–orbit splitting, that is, in MoS₂.

The rest of the paper is organized as follows. In section 2 we introduce the setup and the ingredients of our model: the envelope function approximation, the spin–orbit and Zeeman interactions, the impurity matrix elements, and the electrical drive. In section 3 we describe intervalley scattering between the electronic states, and the important differences between the Mo-type and S-type impurities. In section 4 we introduce the spin-valley qubit. In section 5 we derive the in-plane g-factor of the spin-valley qubit and in section 6 we describe electrically driven Rabi oscillations of the qubit. In section 7 we quantify our generic results in the special case of MoS₂, and conclude in section 8.

The setup we consider is shown in figure 1. The material hosting the spin-valley qubit is a ML-TMDC, e.g. MoS₂, lying in the x-y plane. We consider an electrostatically defined QD, where a single electron is confined with a parabolic, cylindrically symmetric potential

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:vconf} \mathcal{V}_{\rm conf}(x,y)=\frac{1}{2}m^*\omega_0^2\left[(x-x_0)^2+(y-y_0)^2\right], \nonumber \end{align} \tag{ 1 }$$

which is centered at the position $(x_0, y_0)$. Here, m^* is the effective mass of the electron, and $\omega_0$ is the angular frequency of the confinement. The typical length scale of the confinement is the oscillator length, defined as $\newcommand{\e}{{\rm e}} \ell = \sqrt{\hbar / (m^*\omega_0)}$. In principle, the QD can be formed either in the conduction band or in the valence band—throughout the paper we use a terminology corresponding to the conduction band. Here and henceforth, calligraphic font (e.g. ${\mathcal V}_{\rm conf}$) is used to denote real-space Hamiltonians.

In the absence of spin–orbit interaction, magnetic fields, and impurites, the real-space single-electron Hamiltonian of this system reads $\mathcal{H}_{\rm QD} = \mathcal{K} + \mathcal{V}_{\rm cr} + \mathcal{V}_{\rm conf}$, where the first two terms are the kinetic energy and the crystal potential, respectively. The confinement potential $\mathcal{V}_{\rm conf}$ varies slowly in space ($\newcommand{\e}{{\rm e}} \ell \gg a$, where a is the lattice constant), therefore we can build the description of the electronic states upon the envelope-function Hamiltonian

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{envelope} {{\mathcal H}}_{\rm EF}=-\frac{\hbar^2}{2m^*}\left(\frac{\partial^2}{\partial x^2}+\frac{\partial^2}{\partial y^2}\right)+\mathcal{V}_{\rm conf}(x,y). \nonumber \end{align} \tag{ 2 }$$

The eigenfunctions of ${{\mathcal H}}_{\rm EF}$ have the form $\Phi_{nm}(x, y)=\phi_n (x-x_0)\phi_m(y-y_0)$, where $\phi_{n}$ is the nth 1D harmonic-oscillator eigenstate, and the corresponding energy eigenvalue is $(n+m+1)\hbar\omega_0$.

Having the envelope functions $\Phi_{nm}$ at hand, we can express the corresponding real-space wave functions as

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\ket}[1]{|#1\rangle} \displaystyle \label{nmvs_function} \ket{nmvs}=\sum_{\boldsymbol{k}}c_{\boldsymbol k}^{nm}\ket{\Psi^{(1)}_{c,v\boldsymbol{K}+\boldsymbol{k}}} \chi_s. \nonumber \end{align} \tag{ 3 }$$

Here, the spin (valley) degree of freedom is incorporated via the spin quantum number $s\in (\uparrow , \downarrow)$ (valley index $\newcommand{\e}{{\rm e}} v \in (K, K') \equiv (+1, -1)$), $c_{\boldsymbol k}^{nm}$ is the Fourier transform of the envelope function $\Phi_{nm}(x, y)$, and $\chi_s$ is the spin wave function. Furthermore, $\renewcommand{\ket}[1]{|#1\rangle} \ket{\Psi^{(1)}_{c, \boldsymbol{k}}}$ is the conduction-band Bloch function at the $\boldsymbol{k}$ point of the Brillouin zone, approximated up to first order in $k\cdot p$ theory (see details in appendices A and C). Note that $c_{\boldsymbol k}^{nm}$ is a well localized function around $\boldsymbol{k}=0$, so, e.g. the state $\renewcommand{\ket}[1]{|#1\rangle} \ket{nmKs}$ dominantly contains those Bloch functions whose wave vectors are within a small  ∼$\newcommand{\e}{{\rm e}} 1/\ell$ radius of the valley K. The normalization condition is $\sum_{\boldsymbol k} | c_{\boldsymbol k}|^2 =1$. In what follows, we will exploit time reversal symmetry to enforce $\Psi^{(1)}_{c, \boldsymbol K'}(\boldsymbol r) = \left[\Psi_{c, \boldsymbol K}^{(1)}(\boldsymbol r)\right]^*$.

Beside the Hamiltonian $\mathcal{H}_{\rm QD}$, which incorporates the kinetic energy, crystal potential and confinement potential discussed above, we now add to our model the spin–orbit interaction, the magnetic field, the electric field, and the impurity. To formulate this complete Hamiltonian, we use the basis defined in equation (3). In this basis, the real-space Hamiltonian $\mathcal{H}_{\rm QD}$ is represented as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{HamiltonianQD} H_{\rm QD}=\hbar\omega_0\left(a_x^{\dagger}a_x+a_y^{\dagger}a_y+1\right), \nonumber \end{align} \tag{ 4 }$$

with the creation operators defined via

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\ket}[1]{|#1\rangle} \renewcommand{\bra}[1]{\langle#1|} \displaystyle \label{eq5} a_x^{\dagger} = \sum_{nmvs}\sqrt{n+1}\ket{(n+1)mvs}\bra{nmvs}, \nonumber \end{align} \tag{ 5 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\ket}[1]{|#1\rangle} \renewcommand{\bra}[1]{\langle#1|} \displaystyle a_y^{\dagger} = \sum_{nmvs}\sqrt{m+1}\ket{n(m+1)vs}\bra{nmvs}. \nonumber \end{align} \tag{ 6 }$$

Even if the magnetic field is zero, the fourfold degenerate orbital states are split due to the intrinsic spin–orbit interaction, that plays a significant role in ML-TMDC. Typically, spin–orbit splitting $\Delta_{{\rm SO}}$ is expected to exceed (or at least be comparable with) the orbital level spacing $\hbar\omega_0$ in a QD [35]. We describe the spin–orbit splitting using the Hamiltonian

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{so} H_{{\rm SO}}=-\frac{\Delta_{{\rm SO}}}{2}\tau_3 s_z, \nonumber \end{align} \tag{ 7 }$$

where $\tau_3$ is the third Pauli matrix in valley space, defined as $\renewcommand{\bra}[1]{\langle#1|} \renewcommand{\ket}[1]{|#1\rangle} \tau_3 = \sum_{nmv}v\ket{nmv}\bra{nmv}$, and s_x, s_y and s_z are spin Pauli matrices. As a result of the spin–orbit interaction, the fourfold degeneracy of the orbital ground state is split to two Kramers doublets by an energy $\Delta_{\rm SO}$, in a way that the ground-state Kramers doublet is labelled with valley and spin indices as $\renewcommand{\ket}[1]{|#1\rangle} \ket{00K\uparrow}$ and $\renewcommand{\ket}[1]{|#1\rangle} \ket{00K'\downarrow}$. The electron can occupy some superposition within the two-dimensional subspace of the lower-energy Kramers doublet, i.e. the electron represents a Kramers qubit or spin-valley qubit.

The external magnetic field ${\boldsymbol B} =$$(B_\perp\cos{\varphi_B}, B_\perp\sin{\varphi_B}, B_z)$ is represented by the following Hamiltonian

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{magnetic} H_B = \frac 1 2 \mu_B g_s \boldsymbol{B}\cdot\boldsymbol{s} + \frac 1 2 \mu_B g_v B_z \tau_3, \nonumber \end{align} \tag{ 8 }$$

where the terms describe the spin and valley Zeeman splitting. Furthermore, $g_s\approx 2$ is the spin g-factor, $g_v$ is the material dependent valley g-factor, and $\mu_B$ is the Bohr-magneton. In the presence of B_z, the spin-valley qubit basis states (i.e. the lower-energy Kramers doublet $\renewcommand{\ket}[1]{|#1\rangle} \ket{00K\uparrow}$ and $\renewcommand{\ket}[1]{|#1\rangle} \ket{00K'\downarrow}$) are split by the energy $(g_v+g_s)\mu_B B_z$, so the corresponding out-of-plane g-factor $g_{zz}=g_v+g_s$ is finite and also material dependent. On the other hand, in an in-plane magnetic field, the spin-valley qubit remains degenerate, that is, the in-plane g-factor ($g_{xx}=g_{yy}$) in a clean ML-TMDC QD is zero. The simple explanation behind this is that the Kramers-doublet have different spin and valley quantum numbers, and the in-plane B-field does not couple the different valleys.

Note that in equation (8), we neglect orbital effects of the magnetic field beyond the valley Zeeman effect, as those would not contribute to our results within the third-order perturbative description we will apply below. Note, however, that the omitted orbital effects might be important in slightly different settings, e.g. when the spin-valley physics of excited QD orbitals is described.

The impurity is modelled as a spin-independent scattering centre, described by an electrostatic potential $\mathcal{U}_{\rm imp}(\boldsymbol r)$. We choose our reference frame such that the impurity is located at the origin; therefore, the relative position of the impurity and the QD centre is represented by the QD centre location $\boldsymbol r_0 = (x_0, y_0)$ defined above. Because of the short-range character, the impurity couples different orbitals and valley effectively, hence it is expressed in our basis of equation (3) as a dense matrix:

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\ket}[1]{|#1\rangle} \renewcommand{\bra}[1]{\langle#1|} \displaystyle H_{\rm imp} = \sum_{nmv} \sum_{n'm'v'} \sum_s \tilde{\Delta}_{vv'}^{nmn'm'} \ket{nmvs} \bra{n'm'v's}, \nonumber \end{align} \tag{ 9 }$$

where we introduced the QD impurity matrix elements

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\braket}[1]{\langle#1\rangle} \renewcommand{\bra}[1]{\langle#1|} \displaystyle \label{eq:Delta} \tilde{\Delta}_{vv'}^{nmn'm'} = \braket{nmvs | \mathcal{U}_{\rm imp}|n'm'v's}. \nonumber \end{align} \tag{ 10 }$$

In what follows, we will denote the absolute value and the complex phase of the intervalley impurity matrix element as $\Delta_{KK'}^{nmn'm'} > 0$ and $\varphi_{KK'}^{nmn'm'} \in (-\pi, \pi]$.

Finally, to control the spin-valley qubit, we apply an oscillating electric field along the x axis with amplitude E_ac and driving angular frequency ω, which is represented by the real-space Hamiltonian $\mathcal{H}_E = |{\rm e}|xE_{\rm ac}\cos(\omega t)$. In our model, we take this into account via

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{electric} H_E = |{\rm e}|\frac{\ell}{\sqrt{2}} (a_x+a_x^{\dagger}) E_{\rm ac}\cos(\omega t), \nonumber \end{align} \tag{ 11 }$$

where $\rm e$ is the electron charge. This completes our model Hamiltonian $H = H_{\rm QD} + H_{\rm SO} + H_B + H_{\rm imp} + H_E$.

Our goal is to describe the impurity-induced effects (in-plane g-factor and electrically driven qubit Rabi oscillations) using our QD model introduced above. To proceed toward this goal, we need to relate the impurity-induced matrix elements $\tilde{\Delta}_{vv'}^{nmn'm'}$ to the microscopic character of the impurity. As we reveal below, Mo-type impurities and S-type impurities imply qualitatively different intervalley matrix elements, due to symmetry-forbidden intervalley scattering [36, 41, 42].

Combining equations (10) and (3), we express the intervalley QD impurity matrix elements with the bulk Bloch functions $\Psi_{\boldsymbol k}^{(1)}$ as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{KK-general} \tilde{\Delta}_{KK'}^{nmn'm'}= \sum_{\boldsymbol{k},\boldsymbol{k'}} \left(c_{\boldsymbol k}^{nm}\right)^* c_{\boldsymbol k'}^{n'm'} M_{KK'}(\boldsymbol{k},\boldsymbol{k}'), \nonumber \end{align} \tag{ 12 }$$

where we introduced the bulk intervalley impurity matrix elements

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\braket}[1]{\langle#1\rangle} \renewcommand{\bra}[1]{\langle#1|} \displaystyle \label{Mkk} M_{KK'}(\boldsymbol{k},\boldsymbol{k}')=\braket{\Psi^{(1)}_{\boldsymbol{K}+\boldsymbol{k}}|\mathcal{U}_{\rm imp}| \Psi^{(1)}_{\boldsymbol{K'}+\boldsymbol{k}'}}. \nonumber \end{align} \tag{ 13 }$$

These bulk matrix elements and their dependence of the symmetry of the impurity have been characterized, and also quantified for Mo and S vacancies in MoS₂, using density functional theory [41, 42]. Note that it is experimentally established that such defects are dominant in ML-TMDC samples[43–46]. The key findings of the symmetry analysis in [41, 42], and their consequences for our problem, are as follows.

• (i)  
  Mo-type impurity. If the impurity is located at a transition-metal site, e.g. a Mo vacancy in MoS₂, then the bulk intervalley scattering impurity matrix element connecting the K and $K'$ points is finite in general, $M_{KK'}(\boldsymbol 0, \boldsymbol 0) \neq 0$. In this case, since the Fourier transform $c_{\boldsymbol k}^{nm}$ of the envelope function is localized around $\boldsymbol k = 0$, we obtain a finite result if we use the approximation $M_{KK'}(\boldsymbol k, \boldsymbol k') \approx M_{KK'}(\boldsymbol 0, \boldsymbol 0)$ when evaluating equation (12), and this will provide a good approximation. This approximation results in a intervalley QD impurity matrix element

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{KK-M} \tilde{\Delta}_{KK'}^{nmn'm'} = A\Phi^*_{nm}(0,0)\Phi_{n'm'}(0,0)M_{KK'}(\boldsymbol{0}), \nonumber \end{align} \tag{ 14 }$$

  where A is the sample area. Based on the numerical results of [41], for a Mo vacancy in MoS₂ we estimate $M_{KK'}(\boldsymbol 0, \boldsymbol 0) = 145 \, {\rm eV} {\rm \mbox{\AA}}^2/A$.
• (ii)  
  S-type impurity. If the impurity is located at a chalcogen site, e.g. an S vacancy in MoS₂, then the bulk intervalley impurity matrix element, evaluated exactly between K and $K'$, vanishes: $M_{KK'}(\boldsymbol 0, \boldsymbol 0) = 0$. (Note that in the MoS₂ valence band, $M_{KK'}(\boldsymbol{0}, \boldsymbol{0})$ is zero for both Mo and S vacancy.) Numerical results in [41] also reveal that $M_{KK'}(\boldsymbol k, \boldsymbol 0)$ is nonzero, and its absolute value scales linearly with k around $\boldsymbol k = 0$, that is, $|M_{KK'}(\boldsymbol k, \boldsymbol 0)| \propto k$. Using symmetry arguments, we generalize this result in appendix B, where we show that the bulk intervalley matrix element can be described for short wave vectors as

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:sapprox} M_{KK'}(\boldsymbol k,\boldsymbol k') \approx \boldsymbol v \cdot (\boldsymbol k - \boldsymbol k'), \nonumber \end{align} \tag{ 15 }$$

  where $\boldsymbol v = \gamma (1, -i) /A$ and $\gamma > 0$. Based on the numerical results of [41], for an S vacancy in MoS₂ we estimate $\gamma = 15~{\rm eV}~\mathring{\rm A}^3$. The approximation (15) is then used for evaluating the intervalley QD impurity matrix element in equation (12), yielding

  $$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{KK-X} \tilde{\Delta}_{KK'}^{nmn'm'} =\sum_{\boldsymbol{k},\boldsymbol{k'}} \left(c_{\boldsymbol k}^{nm}\right)^* c_{\boldsymbol k'}^{n'm'} \boldsymbol{v}\cdot\left(\boldsymbol{k}-\boldsymbol{k'}\right). \nonumber \end{align} \tag{ 16 }$$

Due to the simple form of the harmonic-oscillator envelope functions, the intervalley QD impurity matrix elements (14) and (16) can be evaluated analytically. Note that here and henceforth we consider the case of a single impurity in the QD. It is straightforward to generalize our results to an impurity ensemble and provide a statistical description, see, e.g. corresponding work in the context of carbon nanotubes [15, 47, 48].

Symmetry-forbidden intervalley scattering, characteristic of S vacancies as described above, can also appear for other short-range impurities. A summary of our findings is shown in table 1. Naturally, for substitutional atoms replacing Mo (S) in the lattice, we expect to find a nonzero (zero) direct intervalley impurity matrix element, as discussed in (i) ((ii)) above, as long as the substitutional atom does not change the corresponding symmetries. Furthermore, we extend the symmetry analysis discussed in this section and in appendix B for the case of adatoms. The four typical high-symmetry locations for adatoms are hollow-site (adatom on the out-of-plane axis piercing the center of a Mo–S hexagon), bridge (adatom in the plane of the S–Mo–S bonds in a given unit cell), atop-Mo (adatom on the out-of-plane axis piercing a Mo atom) and atop-S (adatom on the out-of-plane axis piercing a S atom). In case of an atop-Mo (atop-S) adatom, intervalley scattering is allowed (forbidden), see (i) ((ii)). For a bridge adatom, we find that even though it does preserve a symmetry, namely the S–Mo–S plane within the unit cell remains a mirror plane, this does not give a restriction on the intervalley matrix element. For a hollow-site adatom however, the three-fold rotational symmetry around the axis, together with the symmetry of the conduction-band wave function, implies that direct intervalley scattering is symmetry-forbidden, as in case (ii) above.

We now have all the elements to characterize a spin-valley qubit in a ML-TMDC quantum dot. We consider the perturbative case, when the basis states of the spin-valley qubit are energetically well separated from the other states, i.e. the coupling matrix elements between the qubit and the other states are much smaller than their energy difference:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{hierarchy} \Delta_{{\rm SO}}, \hbar\omega_0 \gg \left\{\begin{array}{@{}c@{}} ~g_s\mu_BB_x \nonumber \\ ~g_v\mu_BB_z \nonumber \\ ~\Delta_{KK'}^{nmn'm'} \nonumber \\ ~|{\rm e}| E_{\rm ac} \ell. \end{array}\right. \nonumber \end{align} \tag{ 17 }$$

Therefore we treat $H_0 = H_{\rm QD} + H_{{\rm SO}}$ as the unperturbed Hamiltonian and $H_1 = H_B+H_{\rm imp}+H_{E}$ as the perturbation. The two lowest-energy eigenstate of the unperturbed Hamiltonian H₀ are $\renewcommand{\ket}[1]{|#1\rangle} \ket{00K\uparrow}$ and $\renewcommand{\ket}[1]{|#1\rangle} \ket{00K' \downarrow}$; we call these the unperturbed spin-valley qubit basis states. Below, we consider the undriven case H_E  =  0 to derive the in-plane g-factor, and the driven case with finite H_E to describe the electrically induced qubit Rabi oscillations. In both cases, we apply Schrieffer–Wolff perturbation theory [49, 50] to derive a $2\times 2$ effective Hamiltonian for the spin-valley qubit. In both cases, the perturbed qubit basis state associated to $\renewcommand{\ket}[1]{|#1\rangle} \ket{00K\uparrow}$ and $\renewcommand{\ket}[1]{|#1\rangle} \ket{00K'\downarrow}$ will be denoted by $\renewcommand{\ket}[1]{|\!\!#1\rangle} \ket{\Uparrow}$ and $\renewcommand{\ket}[1]{|\!\!#1\rangle} \ket{\Downarrow}$, respectively.

In a clean dot, the in-plane magnetic field does not couple the basis states of the spin-valley qubit, so the in-plane g-factor g_xx is zero. However, the impurity-induced intervalley matrix element combined with a finite in-plane magnetic field does couple the qubit basis states via intermediate states at higher energies, see figure 2(a). Therefore, in the presence of impurities, a finite in-plane g-factor is expected for the spin-valley qubit.

Zoom In Zoom Out Reset image size

To evaluate the in-plane g-factor, we use second-order Schrieffer–Wolff perturbation theory [49]. From that, we obtain the effective qubit Hamiltonian

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H_{\rm q} = H_{\rm q}^{(1)} + H_{\rm q}^{(2)}, \label{eq:h2bparallel} \nonumber \end{align} \tag{ 18 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:hq1} H^{(1)}_{\rm q} = \frac{1}{2}\mu_B(g_s+g_v)B_z\sigma_z , \nonumber \end{align} \tag{ 19 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H^{(2)}_{\rm q} = -\frac{\mu_B g_s B_\perp \Delta^{0000}_{KK'}}{\Delta_{\rm SO}} \sigma(\varphi^{0000}_{KK'}+\varphi_B), \label{eq:hq2} \nonumber \end{align} \tag{ 20 }$$

$\sigma(\varphi)=\sigma_x \cos \varphi - \sigma_y \sin \varphi$, and $\sigma_x$, $\sigma_y$, $\sigma_z$ are Pauli matrices acting on the perturbed qubit basis states $\renewcommand{\ket}[1]{|\!\!#1\rangle} \ket{\Uparrow}$ and $\renewcommand{\ket}[1]{|\!\!#1\rangle} \ket{\Downarrow}$. The second-order term $H^{(2)}_{\rm q}$ in equation (20) is interpreted as a Zeeman interaction of the qubit with the in-plane magnetic field $B_\perp$, which is turned on by the presence of the impurity. Figure 2(a) illustrates one term in the perturbative sum that generates $H^{(2)}_{\rm q}$, i.e. one path contributing to this interaction: the lowest two blue lines represent the unperturbed qubit basis states, other blue lines represent other eigenstates of the unperturbed Hamiltonian, whereas the gray arrows represent perturbation matrix elements connecting those states.

From equation (18), we see that the presence of the impurity does not affect the out-of-plane g-factor g_zz, at least in this order of perturbation theory. Furthermore, the in-plane g-factor is expressed as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Gxx} g_{xx} = 2 g_s \frac{\Delta_{KK'}^{0000}}{\Delta_{\rm SO}}. \nonumber \end{align} \tag{ 21 }$$

Note that this is essentially the same result as obtained by Flensberg and Marcus [18] for carbon nanotubes, see their equation (5).

In the rest of this section, we consider impurities that preserve the threefold rotational symmetry of the lattice, either with a rotation axis containing an Mo atom (e.g. a Mo vacancy), or with a rotation axis containing an S atom (e.g. an S vacancy). We show that the dependence of the in-plane g-factor g_xx on the location $\boldsymbol r_0$ of the QD center with respect to the vacancy is qualitatively different for the Mo and S types.

The in-plane g-factor is governed by the intervalley impurity matrix element $\Delta_{KK'}^{0000}$. For a Mo-type impurity, equation (14) implies that this matrix element inherits the spatial dependence of the squared wave function at the impurity position:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{KK_MC} {\Delta_{KK'}^{0000}}=\frac{A}{\ell^2\pi}e^{-\frac{r_0^2}{\ell^2}}M_{KK'}(\boldsymbol{0},{\bf 0}). \nonumber \end{align} \tag{ 22 }$$

Due to equation (21), g_xx inherits the same Gaussian spatial dependence:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:gxxmo} g^{{\rm Mo}}_{xx} = g_{xx,{\rm max}}^{\rm Mo} e^{-\frac{r_0^2}{\ell^2}}, \nonumber \end{align} \tag{ 23 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{gmaxM} g_{xx, {\rm max}}^{\rm Mo}=\frac{2g_sAM_{KK'}(\boldsymbol{0},{\bf 0})}{\ell^2\pi\Delta_{\rm SO}}. \nonumber \end{align} \tag{ 24 }$$

To evaluate $\Delta_{KK'}^{0000}$ for the case of an S-type impurity, we use equation (16), and there we need the Fourier components of the ground-state QD orbital:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{ccoe} c_{\boldsymbol{k}}^{00}=\frac{2\sqrt{\pi}\ell}{\sqrt{A}}e^{-\frac{\ell^2|\boldsymbol k|^2}{2}}e^{-i \boldsymbol{r}_0 \cdot \boldsymbol k}. \nonumber \end{align} \tag{ 25 }$$

(Recall that our reference frame is chosen such that the impurity is located at the origin and the QD center is displaced to $\boldsymbol r_0$.) Then, the intervalley matrix element for the S-type impurity reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{KK_XC} {\Delta_{KK'}^{0000}}= \frac{2}{\pi} \frac{\gamma}{\ell^4}r_0e^{-\frac{r_0^2}{\ell^2}}, \nonumber \end{align} \tag{ 26 }$$

which implies via equation (21) that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle g_{xx}^{\rm S} = g_{xx,{\rm max}}^{\rm S} \sqrt{2 e } \frac{r_0}{\ell} e^{-\frac{r_0^2}{\ell^2}}, \nonumber \end{align} \tag{ 27 }$$

with

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{gmaxX} g_{xx,{\rm max}}^{\rm S}= \frac{2\sqrt{2}}{\pi \sqrt{e}} \frac{g_s\gamma}{\ell^3\Delta_{\rm SO}}. \nonumber \end{align} \tag{ 28 }$$

Clearly, the dependencies of $g_{xx}^{\rm Mo}$ and $g_{xx}^{\rm S}$ on r₀ are qualitatively different, which is illustrated in figures 3(a) and (b). Both plots show circular symmetry: the results depend only on the length of $\boldsymbol r_0$ but not on its direction, see equation (22) and (26). In case of the Mo-type impurity, figure 3(a), the in-plane g-factor decreases monotonically as r₀ is increased, following the Gaussian probability density of the envelope function at the position of the impurity. However, the in-plane g factor due to an S-type impurity, figure 3(b), is different, it is zero if the impurity is at the QD center, and reaches its maximal value $g_{xx, {\rm max}}^{\rm S}$ if the impurity is at a distance $\newcommand{\e}{{\rm e}} r_0 = \ell/\sqrt{2}$ from the QD center.

Zoom In Zoom Out Reset image size

These results imply that in principle, the in-plane g-factor in a MoS₂ QD with a short-range impurity can be tuned electrically, by reshaping or replacing the QD as the voltages on the confinement gates are changed. For example, in the presence of an S-type impurity, the impurity-induced intervalley coupling, and thereby the corresponding decoherence processes, could be suppressed by electrostatically tuning the QD to $(x_0, y_0) = 0$, i.e. without spatially separating the impurity and the electron. Also, the dependence of the in-plane g factor on the impurity position is qualitatively different for Mo-type and S-type impurities. This implies that experiments, where, e.g. magneto-transport spectroscopy maps the g factor as a function of QD shape and location, could reveal information about the type, number, and position of the impurities in the QD.

Here we show that the spin-valley qubit defined in the MoS₂ QD can be coherently controlled by an ac electric field, if a tilted magnetic field is applied, and at least one short-range impurity is present in the dot. In this case, the ac electric field induces qubit Rabi oscillations, when the qubit evolves coherently and cyclically between the two basis states $\renewcommand{\ket}[1]{|\!\!#1\rangle} \ket{\Uparrow}$ and $\renewcommand{\ket}[1]{|\!\!#1\rangle} \ket{\Downarrow}$, if the driving frequency ω is chosen to be resonant with the qubit's Larmor frequency $\omega_q$. Our goal is to calculate the Rabi frequency characterizing these oscillations.

Similarly to the preceding section, here we also map our Hamiltonian $H = H_0 + H_1$, acting in an infinite-dimensional Hilbert space, to a two-dimensional effective qubit Hamiltonian. Note that in this case, H₁ is time dependent due to the presence of the ac electric field in H_E(t). We use time-dependent Schrieffer–Wolff perturbation theory [50], treating H₁ as the perturbation, to obtain the qubit Hamiltonian. We find that the leading-order time-dependent term in the effective qubit Hamiltonian appears in the 3rd order of perturbation theory, and reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle H_{\rm q}^{(3)}(t) = \hbar \Omega_{\rm R} \cos(\omega t) \sigma(\varphi^{0010}_{KK'}+\varphi_B), \label{eq:hqubit3} \nonumber \end{align} \tag{ 29 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:rabigeneric} \hbar \Omega_{\rm R} = \sqrt{2} \frac{(|{\rm e}| E_{ac} \ell) (g_s \mu_B B_\perp) \Delta^{0010}_{KK'} }{\Delta_{\rm SO} \hbar \omega_0}. \nonumber \end{align} \tag{ 30 }$$

This coupling appears due to the interplay of the in-plane magnetic field, the ac electric field and the impurity. Figure 2(b) illustrates one of the three-step paths via virtual intermediate states that contribute to this coupling. Note that the fact that the only intervalley impurity matrix element appearing in equation (29) is $\Delta_{KK'}^{0010}$ is due to the choice that the driving electric field is along the x direction and that the confinement potential is parabolic.

We consider the magnetic-field hierarchy $B_z\gg$$B_\perp \frac{\Delta^{0000}_{KK'}}{\Delta_{\rm SO}}$, for which the qubit Larmor frequency is dominated by the out-of-plane component of the magnetic field. In this case, the qubit Larmor frequency is given by $\omega_q \approx \mu_B(g_s+g_v)B_z$. We want to describe the electrically induced dynamics of this qubit with the initial state being $\renewcommand{\ket}[1]{|\!\!#1\rangle} \ket{\Downarrow}$. Upon resonant or almost-resonant driving, $\omega \approx \omega_q$, and assuming that the driving is weak, the dynamics remains in the two-dimensional space of the spin-valley qubit, and the qubit will show simple Rabi oscillations, characterized by the Rabi frequency $\Omega_{\rm R}$ given by equation (30).

From equation (30), the Rabi frequency is proportional to the impurity matrix element $\Delta^{0010}_{KK'}$, which can be calculated for a Mo-type impurity using equation (14). We find that the dependence of the Rabi frequency as a function of the impurity position is given by

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:rabimo} \Omega_{\rm R}^{\rm Mo} = \Omega_{\rm R,max}^{\rm Mo} \sqrt{\frac{e}{2}}\frac{|x_0|}{\ell} e^{-\frac{r_0^2}{\ell^2}}, \nonumber \end{align} \tag{ 31 }$$

where the maximal value of the Rabi frequency is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{rabimaxM} \Omega_{\rm R,max}^{\rm Mo} = \sqrt\frac{2}{e}\frac{|e| E_{\rm ac}g_s \mu_B B_\perp}{\pi\ell\Delta_{\rm SO} \hbar ^2\omega_0} AM_{KK'}(\boldsymbol{0},\boldsymbol{0}). \nonumber \end{align} \tag{ 32 }$$

The dependence of the Rabi frequency equation (31) on the impurity position is shown in figure 3(c). The key features are: (i) The Rabi frequency is maximized when the position of the impurity is $\newcommand{\e}{{\rm e}} \boldsymbol r_0 =(\pm\frac{\ell}{\sqrt{2}}, 0)$. (ii) The Rabi frequency decreases when the impurity is moved outside from the QD. (iii) The Rabi frequency is exactly zero if the impurity is placed in the QD center, or on the line x₀  =  0.

Feature (iii) can be explained as follows. The external ac electric field along the x axis induces a spatial oscillation of the electron in the x direction, therefore the relative position of the impurity and the QD also oscillates with a small $\newcommand{\e}{{\rm e}} \ll \ell$ amplitude. Since, as seen in figure 3(a), the in-plane g factor g_xx depends on the relative position according to the spatial oscillation of the wave function leads to a temporal oscillation of g_xx, driving a Rabi oscillation between spin-valley qubit basis states. However, in the vicinity of the line x₀  =  0, g_xx depends on x₀ quadratically, therefore the ac electric field does not generate a time-dependent g_xx in first order, and therefore the Rabi frequency is zero.

If the impurity in the QD is S-type, then the Rabi frequency is expressed from equations (30), (16), (25), and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle c_{\boldsymbol k}^{10} = \frac{8\sqrt{\pi}i\ell^2}{\sqrt{A}}k_x e^{-\frac{\ell^2|\boldsymbol k|^2}{2}}e^{-i \boldsymbol{r}_0 \cdot \boldsymbol k}, \nonumber \end{align} \tag{ 33 }$$

as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Omega_{\rm R}^{\rm S} = \Omega_{\rm R,max}^{\rm S} |1-2x_0(x_0+iy_0)/\ell^2| e^{-\frac{r_0^2}{\ell^2}}, \nonumber \end{align} \tag{ 34 }$$

where

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{rabimaxX} \Omega_{\rm R,max}^{\rm S} = \frac{2|e| E_{ac}g_s \mu_B B_\perp}{\pi\ell^2\Delta_{\rm SO} \hbar^2 \omega_0} \gamma. \nonumber \end{align} \tag{ 35 }$$

The spatial dependence of the Rabi frequency is plotted in figure 3(d); the main features are as follows. (i) the Rabi frequency reaches its maximal value if the impurity is in the centre of the QD. (ii) There is no coherent transition if the impurity is placed at the position $\newcommand{\e}{{\rm e}} \boldsymbol r_0 =(\pm\frac{\ell}{\sqrt{2}}, 0)$. (iii) The Rabi frequency has a local maximum if the monovacancy is at the position $\newcommand{\e}{{\rm e}} \boldsymbol r_0 =(\pm\sqrt{\frac 3 2 }\ell, 0)$. We note that the simple explanation of the previous paragraph, applied for the Mo-type impurity and based on the adiabatic modulation of the in-plane g-factor, does not explain features (i) and (iii) in the case of an S-type impurity. We anticipate that in this case, the relation between the g-factor and the qubit dynamics can be characterized by a combination of the 'Zeeman-modulation' and 'iso-Zeeman-modulation' mechanisms [51].

In this work, the in-plane g-factor and the electrically driven Rabi oscillations are described in the perturbative regime, see equation (17), where the spin–orbit splitting $\Delta_{\rm SO}$ is a large energy scale. In this regime, the coupling between the spin-valley qubit basis states decreases as the spin–orbit splitting is increased. Among the semiconductor ML-TMD materials, MoS₂ has the smallest spin–orbit splitting [35]. Therefore, this seems to be the material best suited to observe the effects predicted here. Here we quantify the in-plane g-factor and Rabi frequency for a single-electron spin-valley qubit in a ML MoS₂ QD.

We consider a QD with orbital level spacing $\hbar \omega_0 \approx 0.5$ meV. Using $m^* \approx 0.5~m_{\rm e}$, the oscillator length is $\newcommand{\e}{{\rm e}} \ell \approx 17.5$ nm. According to numerical studies of ML MoS₂, the spin–orbit gap is $\Delta_{\rm SO}=3$ meV and the valley g factor is $g_v\approx 4$ [35]. We assume a magnetic field with in-plane component $B_\perp=1$ T, and a driving electric field with amplitude $E_{\rm ac}=10$ kV m⁻¹. The latter value implies that the amplitude of the electrons spatial oscillations is 6 nm, smaller than the oscillator length $\newcommand{\e}{{\rm e}} \ell$, so the perturbative treatment described above is justified.

If the Mo-type impurity is a Mo vacancy, and it is placed at the QD center, then equation (22) combined with the numerical estimate [41] $M_{KK'}({\bf 0}, {\bf 0}) = 145~{\rm eV}~\mathring{\rm A}^2/A$ gives $\Delta_{KK'}^{0000}=1.5$ meV for the intervalley impurity matrix element. This value is larger than the orbital level spacing, in disagreement with the assumption equation (17). Therefore, our results (23) and (31) for the g-factor and Rabi frequency can be applied for a Mo vacancy in MoS₂ only if the distance of the vacancy from the QD centre is larger than 32 nm; in this case the orbital level spacing exceeds the impurity matrix element, and hence a perturbative treatment of the latter is valid. Of course, our results (23) and (31) can also be applied for other Mo-type defects in MoS₂ or in other ML-TMDC materials, if the corresponding impurity matrix elements are below the orbital level spacing.

If the S-type impurity is a S vacancy, then the maximal value of intervalley impurity matrix element according to equation (26) is 0.8 μeV, complying with the perturbative hierarcy of equation (17). Then, the value of the maximal in-plane g factor equation (28) is $g_{xx, {\rm max}}^{\rm S}=10^{-3}$, and the maximal Rabi-frequency from equation (35) is $\Omega_{{\rm R, max}}^{\rm S}=37$ MHz. Note that in this example, also the intravalley QD impurity matrix element is smaller then the orbital level spacing, and therefore using the harmonic-oscillator envelope functions as the orbital basis is a reasonable approximation. (Based on the numerical results of [41], we estimate $M_{KK}({\bf 0}, {\bf 0}) = 15~ {\rm eV}~\mathring{\rm A}^2/A$, implying $\Delta_{KK}^{0000}=150$ μeV for the parameter set considered here.)

In this work, we have proposed and analysed a way to coherently control a single-electron spin-valley qubit, that is defined in a QD in a ML-TMDC material, e.g. MoS₂. Qubit control is performed in a fashion similar to electrically driven spin resonance: a resonant ac electric field drives Rabi oscillations with the help of an in-plane magnetic field and a short-range impurity; the former flips the spin and the latter flips the valley. For the case of a S-type impurity in the QD, we estimated that electrically driven qubit Rabi frequencies of the order of 10–100 MHz can be achieved. We also revealed and discussed the unconventional dependence of the in-plane g factor and the electrically driven Rabi frequency on the impurity position, which arises for S-type impurities due to symmetry-forbidden intervalley scattering. Our results provide design guidelines for efficient electric control of spin-valley qubits in monolayer MoS₂, which are promising building blocks for quantum information processing experiments, expected to show boosted coherence times, especially in isotopically purified samples.

We acknowledge useful discussions with G Burkard, R Frisenda, K Kaasbjerg, J Pető, and L Tapasztó. GSz and AP acknowledge the financial support of the National Research, Development and Innovation Office of Hungary via the National Quantum Technologies Program NKP-2017-00001 and the OTKA Grants 105149, 108676 and 124723. LC acknowledges funding from the European Union's Seventh Framework Programme (FP7/2007-2013) through the ERC Advanced Grant NOVGRAPHENE (GA No. 290846) and the Comunidad de Madrid through the grant MAD2D-CM, S2013/MIT-3007. AP is supported by the ÚNKP-17-4-III New National Excellence Program of the Ministry of Human Capacities of Hungary.

Here, we review the elements of $k\cdot p$ theory that are relevant for the description of conduction-band electrons in the ML-TMDC materials studied here. These considerations are required to describe intervalley scattering by S-type impurities, e.g. an S vacancy in MoS₂.

The single-electron Hamiltonian of the perfect crystal is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{H0crystal} \mathcal{H}_{\rm cr}= \frac{p_x^2+p_y^2+p_z^2}{2m_{\rm e}} + \mathcal{V}_{\rm cr}(\boldsymbol r,z), \nonumber \end{align} \tag{ A.1 }$$

where the crystal potential is denoted by $\mathcal{V}_{\rm cr}(\boldsymbol r, z)$ with $\boldsymbol r=(x, y)$. The eigenstates of ${\mathcal H}_{\rm cr}$ are Bloch states:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{Bloch} \Psi_{\alpha, \boldsymbol k} (\boldsymbol r,z)= \frac{1}{\sqrt{A}} e^{i\boldsymbol k \boldsymbol r} u_{\alpha,\boldsymbol k} (\boldsymbol r,z), \nonumber \end{align} \tag{ A.2 }$$

where α is the band index, A is the sample area, and the lattice-periodic u functions are normalized as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \int_{-\infty}^{\infty} dz \int_{\rm sample} d^2 \boldsymbol r |u_{\alpha, \boldsymbol k}(\boldsymbol r,z)|^2 = A. \nonumber \end{align} \tag{ A.3 }$$

The corresponding eigenvalues (i.e. the band structure) are denoted as $\newcommand{\e}{{\rm e}} \epsilon_{\alpha, \boldsymbol k}$.

Our quantum dot holds a single conduction electron in a spatially slowly varying confinement potential. In this setup, the electronic wave function can be described via Bloch states from the vicinity of the conduction-band edge, that is, from the vicinity of the K and $K'$ points of the Brillouin zone. For crystal momenta around, e.g. K, we can use first-order perturbation theory to express the conduction-band ($\alpha = c$) Bloch states $\Psi_{c, \boldsymbol K +\boldsymbol k}$ with the band-edge Bloch states $\Psi_{\alpha, \boldsymbol K}$—this procedure is called $k\cdot p$ perturbation theory. In this theory, the so-called Kohn–Luttinger basis functions are used to form a matrix representation of the Hamiltonian:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:kohnluttinger} \phi_{\alpha, \boldsymbol K + \boldsymbol k} (\boldsymbol r,z) &= e^{i\boldsymbol k \boldsymbol r} \Psi_{\alpha, \boldsymbol K}(\boldsymbol r,z) \nonumber \\ &= \frac{1}{\sqrt{A}} e^{i\boldsymbol K \boldsymbol r} e^{i\boldsymbol k \boldsymbol r} u_{\alpha \boldsymbol K}(\boldsymbol r,z). \nonumber \end{align} \tag{ A.4 }$$

In the Kohn–Luttinger basis, the crystal Hamiltonian reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{KL_Hamilton} \left(\mathcal{H}_{\rm cr}\right)_{\alpha, \boldsymbol K+\boldsymbol k, \alpha', \boldsymbol K+ \boldsymbol k'} = \epsilon_{\alpha} \delta_{\alpha\alpha'} \delta_{\boldsymbol k \boldsymbol k'} + \delta_{\boldsymbol k \boldsymbol k'} \boldsymbol k \cdot \boldsymbol P_{\alpha\alpha'}, \nonumber \end{align} \tag{ A.5 }$$

where $\newcommand{\e}{{\rm e}} \epsilon_\alpha = \epsilon_{\alpha, \boldsymbol K}$, the fully diagonal kinetic energy term $\hbar^2 k^2 / (2m_{\rm e})$ was omitted, and

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\braket}[1]{\langle#1\rangle} \renewcommand{\bra}[1]{\langle#1|} \displaystyle \label{Pkdotp} \boldsymbol P_{\alpha\alpha'} = \frac 1 A \frac{\hbar}{m_{\rm e}} \braket{u_{\alpha\boldsymbol K}| \boldsymbol p | u_{\alpha'\boldsymbol K}}(1-\delta_{\alpha,\alpha'}) , \nonumber \end{align} \tag{ A.6 }$$

with $\boldsymbol p = (\,p_x, p_y)$.

The next step in $k\cdot p$ theory is to reduce the multi-band problem defined by the crystal Hamiltonian $\mathcal{H}_{\rm cr}$ to the conduction band. For this, second-order Schrieffer–Wolff perturbation theory [49] can be applied, where the fully diagonal first term of equation (A.5) is the unperturbed Hamiltonian, and the second ($k\cdot p$) term of equation (A.5) forms the perturbation. The Schrieffer–Wolff transformation has the form $\tilde {\mathcal{H}}_{\rm cr} = e^{-S} \mathcal{H}_{\rm cr} e^{S}$, where the expression for the anti-Hermitian first-order transformation matrix is [49]

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:schriefferwolff} S_{c\boldsymbol K + \boldsymbol k, \alpha \boldsymbol K + \boldsymbol k'} \approx - \delta_{\boldsymbol k \boldsymbol k'} \frac{\boldsymbol k \cdot \boldsymbol P_{c\alpha}}{\epsilon_{c}-\epsilon_{\alpha}}, \nonumber \end{align} \tag{ A.7 }$$

and the resulting effective conduction-band Hamiltonian reads

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \left(\tilde{\mathcal{H}}_{\rm cr}\right)_{c\boldsymbol K + \boldsymbol k,c\boldsymbol K + \boldsymbol k'} = \boldsymbol k \frac{\hbar^2}{2m^*} \boldsymbol k \, \delta_{\boldsymbol k,\boldsymbol k'}. \label{eq:hkp} \nonumber \end{align} \tag{ A.8 }$$

Here, we introduced the inverse effective mass tensor via

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac 1 {m^*} = \frac{2}{\hbar^2} \sum_{\alpha \neq c} \frac{\boldsymbol P_{c \alpha} \circ \boldsymbol P_{\alpha c} }{\epsilon_{c} - \epsilon_{\alpha}} \nonumber \end{align} \tag{ A.9 }$$

where ° is the dyadic product. Note that since a threefold rotation around an out-of-plane axis is a symmetry of the lattice, the effective mass is isotropic, i.e. this tensor is diagonal.

The effective conduction-band Hamiltonian matrix $\tilde{{\mathcal H}}_{\rm cr}$ is diagonal, so its eigenvectors $c^{(\boldsymbol k)}$ can be labelled by the wave vector $\boldsymbol k$, and have the trivial structure $c^{(\boldsymbol k)}_{\boldsymbol k'} = \delta_{\boldsymbol k, \boldsymbol k'}$. The corresponding eigenvalue is $\hbar k^2/(2m^*)$. To obtain the corresponding Bloch state $\Psi_{c, \boldsymbol K + \boldsymbol k}(\boldsymbol r, z)$, we have to invert the Schrieffer–Wolff transformation,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:swinverse} \Psi_{c,\boldsymbol K + \boldsymbol k} = \sum_{\alpha \boldsymbol k'} \left[ \sum_{\boldsymbol k'' } \left(e^S \right)_{\alpha,\boldsymbol K + \boldsymbol k',c \boldsymbol K + \boldsymbol k''} c^{(\boldsymbol k)}_{\boldsymbol k''} \right] \phi_{\alpha,\boldsymbol K + \boldsymbol k'}. \nonumber \end{align} \tag{ A.10 }$$

Making a linear expansion of the exponential, and using the the first-order result (A.7) for the transformation matrix, we obtain the perturbative first-order approximation $\Psi^{(1)}_{c, \boldsymbol K + \boldsymbol k}$ for the Bloch state $\Psi_{c, \boldsymbol K + \boldsymbol k}$ as

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \Psi^{(1)}_{c,\boldsymbol K + \boldsymbol k} = \phi_{c,\boldsymbol K + \boldsymbol k}+ \sum_{\alpha \neq c} \frac{\boldsymbol k \cdot \boldsymbol P_{\alpha c}}{\epsilon_{c } - \epsilon_{\alpha}} \phi_{\alpha, \boldsymbol K + \boldsymbol k}, \label{eq:perturbativebloch} \nonumber \end{align} \tag{ A.11 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle = e^{i\boldsymbol k \boldsymbol r} \left(\Psi_{c,\boldsymbol K}+ \sum_{\alpha \neq c} \frac{\boldsymbol k \cdot \boldsymbol P_{\alpha c}}{\epsilon_{c} - \epsilon_{\alpha }} \Psi_{\alpha, \boldsymbol K} \right), \label{eq:kpresult} \nonumber \end{align} \tag{ A.12 }$$

where the arguments $\boldsymbol r, z$ are suppressed. This result is used in this work to evaluate the linear-in-momentum component of the intervalley matrix elements of a short-range impurity. Note that this treatment, i.e. taking into account the $k\cdot p$-induced hybridisation of the conduction band with remote bands, is required only for S-type impurities, where a linear-in-$\boldsymbol k$ expansion of the intervalley matrix element is needed to obtain a nonzero result.

For an S-type impurity, the symmetries of the conduction-band-edge Bloch function and the atomic structure pose the constraint (15) on the linear-in-$\boldsymbol k$ contribution of the bulk intervalley impurity matrix element $M_{KK'}(\boldsymbol k, \boldsymbol k')$, defined in equation (13). Here we derive this constraint, with the specification that the Bloch states in the definition (13) are given by the $k\cdot p$ result (A.12).

The generic expression for the linear-in-$\boldsymbol k$ bulk intervalley impurity matrix element, without any symmetry restrictions, is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:mkkgeneral} M_{KK'}(\boldsymbol k,\boldsymbol k') = z_0 + z_x k_x + z_y k_y + z'_x k'_x + z'_y k'_y, \nonumber \end{align} \tag{ B.1 }$$

where the 5 coefficients are complex numbers. It is natural to expect that the rotational symmetry (C₃) provides a relation between $M_{KK'}(C_3 \boldsymbol k, C_3\boldsymbol k')$ and $M_{KK'}(\boldsymbol k, \boldsymbol k')$, and that time-reversal symmetry provides a relation between $M_{KK'}(-\boldsymbol k', -\boldsymbol k)$ and $M_{KK'}(\boldsymbol k, \boldsymbol k')$. We derive these relations here, transform them to relations between the 5 z coefficients, and thereby simplify the general expression in equation (B.1).

To exploit the rotational symmetry C₃, we start from

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\braket}[1]{\langle#1\rangle} \renewcommand{\bra}[1]{\langle#1|} \displaystyle \label{eq:mkkc3} M_{KK'}(C_3 \boldsymbol k, C_3\boldsymbol k') = \braket{\Psi^{(1)}_{c,\boldsymbol K+C_3 \boldsymbol k} | \mathcal{U}_{\rm imp} | \Psi^{(1)}_{c,\boldsymbol K'+C_3 \boldsymbol k'} }. \nonumber \end{align} \tag{ B.2 }$$

On the one hand, from equation (A.12), we know

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:rotation1} \Psi^{(1)}_{c,\boldsymbol K + C_3 \boldsymbol k} = e^{i (C_3 \boldsymbol k) \boldsymbol r} \left[ \Psi_{c,\boldsymbol K}+ \sum_{\alpha \neq c} \frac{(C_3 \boldsymbol k) \cdot \boldsymbol P_{\alpha c}}{\epsilon_{c \boldsymbol K} - \epsilon_{\alpha \boldsymbol K}} \Psi_{\alpha, \boldsymbol K} \right],\nonumber \\ \nonumber \end{align} \tag{ B.3 }$$

and that this is a non-degenerate Bloch state with energy eigenvalue $\newcommand{\e}{{\rm e}} \epsilon_{c, \boldsymbol K + \boldsymbol k}$. On the other hand, we know [36, 41] how the band-edge Bloch states transform under rotation around an S site: $C_3 \Psi_{\alpha, \boldsymbol K} = \omega^{n_\alpha} \Psi_{\alpha, \boldsymbol K}$ with $\omega = e^{i2\pi/3}$ and $n_\alpha \in \{-1, 0, 1\}$, and n_c  =  1 in particular. Combining these with equation (A.12), we find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:rotation2} C_3 \Psi^{(1)}_{c,\boldsymbol K+ \boldsymbol k} = e^{i (C_3 \boldsymbol k) \boldsymbol r} \left[ \omega \Psi_{c,\boldsymbol K}+ \sum_{\alpha \neq c} \frac{\boldsymbol k \cdot \boldsymbol P_{\alpha c}}{\epsilon_{c \boldsymbol K} - \epsilon_{\alpha \boldsymbol K}} \omega^{n_\alpha}\Psi_{\alpha, \boldsymbol K} \right]. \nonumber \end{align} \tag{ B.4 }$$

Importantly, this is also a Bloch state with energy $\newcommand{\e}{{\rm e}} \epsilon_{c, \boldsymbol K + \boldsymbol k}$. Since the Bloch states (B.3) and (B.4) have the same energy and momentum, they must be identical, up to a complex phase factor. This requirement implies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:omegaK} \Psi^{(1)}_{c,\boldsymbol K + C_3 \boldsymbol k} = \omega^{-1} C_3 \Psi^{(1)}_{c,\boldsymbol K + \boldsymbol k}. \nonumber \end{align} \tag{ B.5 }$$

We describe the $K'$-valley Bloch states as time-reversed K-valley Bloch states, $\Psi_{\alpha, \boldsymbol K'} = \Psi_{\alpha, \boldsymbol K}^*$, yielding $C_3 \Psi_{\alpha, \boldsymbol K'} = \omega^{-n_\alpha} \Psi_{\alpha, \boldsymbol K'}$, and

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:omegaKprime} \Psi^{(1)}_{c,\boldsymbol K' + C_3 \boldsymbol k'} = \omega C_3 \Psi^{(1)}_{c,\boldsymbol K' + \boldsymbol k'}. \nonumber \end{align} \tag{ B.6 }$$

Combining equations (B.5) and (B.6) in equation (B.2), and exploiting the rotational symmetry $C_3 \mathcal{U}_{\rm imp} C_3^{-1} = \mathcal{U}_{\rm imp}$ of the impurity potential, results in

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:mkkrot} M_{KK'} (C_3 \boldsymbol k, C_3\boldsymbol k') = \omega^{-1} M_{KK'} (\boldsymbol k, \boldsymbol k'). \nonumber \end{align} \tag{ B.7 }$$

To exploit time reversal symmetry T, which is the complex conjugation in our case, fulfilling T  =  T⁻¹, we start from

$$\begin{align} \newcommand{\e}{{\rm e}} \renewcommand{\braket}[1]{\langle#1\rangle} \renewcommand{\bra}[1]{\langle#1|} \displaystyle \label{eq:mkktrs} M_{KK'}(-\boldsymbol k',-\boldsymbol k) = \braket{\Psi^{(1)}_{c,\boldsymbol K - \boldsymbol k'} | \mathcal{U}_{\rm imp} | \Psi^{(1)}_{c,\boldsymbol K' - \boldsymbol k} }. \nonumber \end{align} \tag{ B.8 }$$

Again, we describe the $K'$-valley Bloch states as time-reversed K-valley Bloch states: $\Psi^{(1)}_{c, \boldsymbol K'-\boldsymbol k} = \left[\Psi^{(1)}_{c, \boldsymbol K + \boldsymbol k}\right]^*$. Using this together with the time-reversal symmetry $T\mathcal{U}_{\rm imp}T^{-1} = \mathcal{U}_{\rm imp}$ of the impurity potential, i.e. that the latter is real-valued, we find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_{KK'}(-\boldsymbol k',-\boldsymbol k) = M_{KK'} (\boldsymbol k, \boldsymbol k'). \nonumber \end{align} \tag{ B.9 }$$

The symmetry-based equations (B.7) and (B.8) must hold for any values of the wave-vector components k_x, k_y, $k'_x$, $k'_y$. Therefore, these equations form a set of 10 complex linear equations for the 5 complex unknown z coefficients. We solve this linear set, and apply the resulting relations in equation (B.1), yielding

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle M_{KK'} (\boldsymbol k, \boldsymbol k') = (k_x - i k_y - k'_x + i k'_y) z_x. \nonumber \end{align} \tag{ B.10 }$$

Clearly, z_x inherits a global complex phase factor from the band-edge Bloch state $\Psi_{c, \boldsymbol K}$, and hence with appropriate choice of global complex phase for $\Psi_{c, \boldsymbol K}$ we can make z_x positive. Hence, we have proven that the form of the bulk intervalley impurity matrix elements follow equation (15), with the identification $z_x = \gamma/A$.

Note that the key difference between the Mo-type and S-type impurities is the transformation factor n_c. I.e. both types of impurities preserve the C₃ rotation around the impurity site as a symmetry of the structure. However, the rotation around the Mo-type impurity transforms [41] the conduction-band-edge state $\Psi^{(1)}_{c, \boldsymbol K}$ as $C_3 \Psi^{(1)}_{c, \boldsymbol K} = \Psi^{(1)}_{c, \boldsymbol K}$, that is, with n_c  =  0, implying that $M_{KK'}({\bf 0}, {\bf 0})$ is finite in general for a Mo-type impurity, and therefore can be used to approximate the intervalley impurity matrix elements between the QD wave functions.

We have described the wave function of the electron in our QD via equation (3), i.e. as a wave packet, formed by the conduction-band Bloch states $\Psi^{(1)}_{c, \boldsymbol k}$ expressed from $k\cdot p$ theory. Here we summarize how equation (3) is obtained from the envelope-function approximation (EFA).

The Hamiltonian $\mathcal{H} = \mathcal{H}_{\rm cr}+ \mathcal{V}_{\rm conf}$ of the electron in the QD is the sum of the crystal Hamiltonian $\mathcal{H}_{\rm cr}$ and the QD confinement potential $\mathcal{V}_{\rm conf}$. The EFA is an extension of $k\cdot p$ theory: it also relies on the Kohn–Luttinger matrix representation of the Hamiltonian, which is now given by equation (A.5) plus $\left(\mathcal{V}_{\rm conf}\right)_{\alpha\boldsymbol k, \alpha' \boldsymbol k'}$. Exploiting the spatially slowly varying nature of $\mathcal{V}_{\rm conf}$ and the orthogonality of the $u_{\alpha \boldsymbol K}$ functions of different bands, the interband matrix elements are neglected, $\left(\mathcal{V}_{\rm conf}\right)_{\alpha\boldsymbol k, \alpha' \boldsymbol k'} = \delta_{\alpha\alpha'}\left(\mathcal{V}_{\rm conf}\right){\hspace{0pt}}^{(\alpha)}_{\boldsymbol k \boldsymbol k'}$. Furthermore, intervalley matrix elements of $\mathcal{V}_{\rm conf}$ are neglected. This approximation is justified by the numerical results of [36], which imply that in our example, the order of magnitude of the confinement-induced intervalley matrix element does exceed 1 neV, it remains well below the impurity-induced intervalley matrix element of a single S vacancy (∼1 μeV, see our section 7).

Similarly to $k\cdot p$ theory, the next step in the EFA is to reduce the multi-band problem of $\mathcal{H}$ to the conduction band via the Schrieffer–Wolff transformation, treating the confinement potential as part of the perturbation. Importantly, since the interband matrix elements of the confinement potential are neglected, the interband perturbation remains the same as in the absence of the confinement potential, which implies that the Schrieffer–Wolff transformation matrix S has exactly the same form as in $k\cdot p$ theory, i.e. equation (A.7). The transformation results in the following effective conduction-band Hamiltonian (see equation (A.8)):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \tilde{\mathcal{H}}_{c\boldsymbol K + \boldsymbol k,c\boldsymbol K + \boldsymbol k'} = \frac{\hbar^2 k^2}{2m^*} \, \delta_{\boldsymbol k,\boldsymbol k'} + \left(\mathcal{V}_{\rm conf} \right)^{(c)}_{\boldsymbol K + \boldsymbol k,\boldsymbol K+ \boldsymbol k'}. \label{eq:hefkspace} \nonumber \end{align} \tag{ C.1 }$$

Recall that the effective mass tensor is isotropic; this fact was used in equations (1) and (2).

Solutions $(c, E)$ of the matrix Schrödinger equation

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sum_{\boldsymbol k'} \tilde{\mathcal{H}}_{c\boldsymbol K + \boldsymbol k,c\boldsymbol K + \boldsymbol k'} c_{\boldsymbol k'} = E c_{\boldsymbol k}, \nonumber \end{align} \tag{ C.2 }$$

are usually derived by Fourier transforming this equation to a real-space envelope-function Schrödinger equation, yielding equation (2), finding the solution $(\Phi, E)$ of the latter, and Fourier transforming the envelope function $\Phi(x, y)$ to momentum space, yielding the coefficients $c_{\boldsymbol k}$.

To express the real-space wave function $\psi(\boldsymbol r, z)$ corresponding to the eigenvector c of the transformed Hamiltonian $\tilde{\mathcal{H}}$, we have to invert the Schrieffer–Wolff transformation (see equation (A.10)):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \psi = \sum_{\alpha,\boldsymbol k} \left[ \sum_{\boldsymbol k' } \left(e^S \right)_{\alpha,\boldsymbol K + \boldsymbol k,c \boldsymbol K + \boldsymbol k'} c_{\boldsymbol k'} \right] \phi_{\alpha,\boldsymbol K + \boldsymbol k}. \nonumber \end{align} \tag{ C.3 }$$

Using a linear expansion of the Schrieffer–Wolff transformation matrix in the perturbation, we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \psi = \sum_{\boldsymbol k} c_{\boldsymbol k} \Psi^{(1)}_{c,\boldsymbol K + \boldsymbol k}, \label{eq:realspacewfn} \nonumber \end{align} \tag{ C.4 }$$

which is then utilized in the main text as equation (3). In words, equation (C.4) expresses the fact that each electronic real-space energy eigenfunction of the QD is a packet of Bloch waves from the vicinity of the K point, and the coefficients $c_{\boldsymbol k}$ of this wave packet can be obtained from the envelope-function Schrödinger equation.

In the main text, we have presented the effective spin-valley qubit Hamiltonian that is required to describe the in-plane g factor and the electrically driven Rabi oscillations. For completeness, here we present the third-order effective Hamiltonian, including the third-order time-independent term that was not discussed in the main text.

Recall that performing a second-order (time-independent) Schrieffer–Wolff transformation yields $H_{\rm q} = H^{(1)}_{\rm q} + H^{(2)}_{\rm q}$ as shown by equations (18), (19), and (20). Performing a third-order time-dependent Schrieffer–Wolff transformation [50] we find

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:tdswresult} H_{\rm q} = H_{\rm q}^{(1)} + H^{(2)}_{\rm q} + H^{(3)}_{\rm q} + H^{(3)}_{\rm q}(t), \nonumber \end{align} \tag{ D.1 }$$

where the first- and second-order terms are the same as above, the time-dependent third-order term is given in equation (29), and the time-independent third-order term is

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eq:h3e} H_{\rm q}^{(3)} =& \frac 1 4 \mu_B^3 g_s^3 \frac{B_\perp^2 B_z}{\Delta_{{\rm SO}}^2} \sigma_z \nonumber \\ & - g_v \mu_B B_z\sum_{nm}\frac{\left(\Delta_{KK'}^{00nm}\right)^2}{(\Delta_{{\rm SO}}+(n+m)\hbar \omega_0)^2}\sigma_z\nonumber \\ &+2 \frac{g_s \mu_B B_\perp }{\Delta_{{\rm SO}} \hbar \omega_0} \sum_{{nm \atop (nm)\neq(00)}}\frac{1}{n+m} \left[ \Delta^{00nm}_{KK'} \Delta^{00nm}_{KK} \right] \nonumber \\ & \times\sigma(\varphi_{KK'}^{00nm}+\varphi_B). \nonumber \end{align} \tag{ D.2 }$$

As in the main text, $\sigma_x$, $\sigma_y$, $\sigma_z$ are Pauli matrices acting on the perturbed states $\renewcommand{\ket}[1]{|\!\!#1\rangle} \ket{\Uparrow}$ and $\renewcommand{\ket}[1]{|\!\!#1\rangle} \ket{\Downarrow}$, see section 4. The three subsequent terms describe the nonlinear Zeeman effect, and a higher-order correction of the out-of-plane and in-plane g-factors, respectively.